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Abstract 

We present results of Diffusion Monte Carlo calculations for a system of solid ortho-D2 at dif- 
ferent densities, for pressure ranging from up to 350MPa. We compare the equation of state 
obtained using two of the most used effective intermolecular potentials, i. e. the Silvera-Goldman 
and the Buck potentials, with experimental data, in order to assess the validity of the model inter- 
actions. The Silvera-Goldman potential has been found to provide a satisfactory agreement with 
experimental results, showing that, as opposed to what recently found for P-H2, three-body forces 
can be efficiently accounted for by an effective two-body term. 



PACS numbers: 67.80.Mg,67.80.-s 



INTRODUCTION 



Solid molecular hydrogen has been the subject of intense experimental and theoretical 
study over the years. However, the properties of bulk isotopic equivalents (like HD or D 2 ) 
were not investigated in comparable detail. 

It is a common practice to describe molecules of para-hydrogen and isotopic equivalents 
as point particles interacting via a model potential which in some case includes terms which 
effectively account for three-body interactions, which are expected to be not negligible in 
these systems. The popular Silvera-Goldman (SG)|l| potential, which includes a 1/r 9 term 
which mimicks a triple-dipole interaction, and the Buck j^, [3] potential, which is instead 
a pure two-body interaction, have been employed indifferently in theoretical works about 
P-H2, while in the case of 0-D2 only the first one has been commonly used. However, in the 
case of p-H 2 , we recently discovered by means of accurate Quantum Monte Carlo simulations 
that an explicit three-body term is necessary in order to correctly reproduce experimental 
data on pressure Q], . It is therefore interesting to investigate if the different mass plays a 
role in determining the necessity of including an explicit three-body force in the description 
of o-D 2 . As already mentioned, both the Buck and the SG interactions are isotropic and 
therefore either completely neglect, or treat in an effective way the possible three body 
interactions acting among the molecules. This commonly used approximation is in principle 
justified by the fact that at T=0K the molecules of hydrogen and deuterium are considered 
to be in their rotational ground-state (J=0) and therefore the average interaction must have 
spherical symmetry. Three body interactions should therefore be a second order effect. 

In this work we present the results of a systematic study of the equation of state of solid 
0-D2 for pressures ranging from (and slightly below) up to 350MPa performed by means 
of the Diffusion Monte Carlo algorithm, which gives for many Bosons the exact eigenvalue 
of the system, depending only on the chosen Hamiltonian. We compared the results of our 
simulations with the experimental measurements of the pressure. Only the SG potential 
is found to be capable to give a satisfactory description of the pressure-volume curve in 
the studied range of density. We can safely attribute this result to the presence in the 
potential of the repulsive extra-term which treats in an effective way the possible three body 
interactions acting among the molecules. 

The next section of the paper will describe the methods used in our analysis. Section III 
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will present the results on energies, pressures and and other ground state properties in o-D 2 . 
Section IV is devoted to conclusions. 



METHODS 

We model a system of solid o-D 2 as a set of N point particles described by the following 
Hamiltonian: 

« = 4^ + E^(4 (!) 

1=1 Kj 

where R = {ri . . . tn} stands for the 3N coordinates of the molecules and v 2 (R) is the 
intermolecular potential. As already mentioned, we use two different molecule-molecule 
interaction t> 2 , i.e. the Silvera-Goldman and the Buck potentials. 

Diffusion Monte Carlo is a stochastic algorithm projecting the ground state of a many- 
body Bosonic system by starting from a variational ansatz for the many-body wavefunction 
^t(R) and evolving the function f(R, r) = ^/t(R)^(R, t) in imaginary time t — it accord- 
ing to the time-dependent Schrodinger equation: 

df(R,r) 1 2 



dr 2 



--V7(^,r) + V- F(R)f(R,r) 



+ 



E L (R)-E T f(R,r). (2) 



Here the quantum force F is defined as: 

F (fl) = vM» T (*)] = T^ (3) 

and the coefficient of the branching term is defined in terms of the local energy: 
Under the assumption that 

(*(i?,O)|0 o ) = (* T (i2)|0o>^O, (5) 

in the limit of r — > +oo the distribution f(R,r) becomes proportional to ^t{R)4>o{R)- 

In the actual simulation one evolves a population of N w points in configuration space 
(walkers) drawn from the initial distribution f(R,0) = \^t\ 2 according the approximate 



3 



importance sampled Green's function (for details about the sampling algorithm see e. g. 
Ref. E): 

G(R, R', At) = G (R, R', A ^)|^y (6) 

where 

G (R, R', Arje-T^Me^e-T^ 11 '). (7) 

The propagation has to be iterated in order to achieve a propagation for an imaginary time 
t such that the walkers evolved in the Monte Carlo are asymptotically sampled from ^t4>o 
In this situation the expectation value of the Hamiltonian itself computed as the average of 
the local energy riip{R) /ip{R) on the sampled walkers becomes: 



J dR^ T (R)MR)^T(R)/^T(R) „ 

{"-) — r j n,T, TTvTi 77v\ — -^0, 



(8) 



fdR* T (R)MR) 

the lowest eigenvalue of the Hamiltonian. The computed estimate is exact at order At and 
the results are are largely independent of the choice of the trial wave function. Nevertheless, 
^>t has be variationally optimized, in order to reduce the variance on the estimate of the 
eigenvalue in the DMC calculation. We chosen the popular Jastrow-Nosanow form 



M I N f 

V>r(-R) = n ex P 7; u ( r ij) II ex P [ ~ Cfa - Si) 

i<j -I i=l 



(9) 



where the {Si} are the coordinates of the lattice sites around which the molecules are 
confined. 

The possibility of neglecting the symmetrization under exchange of molecules in the trial 
wavefunction has been already discussed in the case of solid p-H 2 [^], justifing this ansatz 
also for solid 0-D2 which is expected to be even more localized than and P-H2. 

A series expansion of the Jastrow factor in equation was used in order to improve the 
description of the two-body correlations in the system. The function u(r) is the sum of a 
standard McMillan pseudopotential and of a correction term which has been expanded in 
terms of suitable basis functions: 



u{r) 



(10) 



The basis functions \n (about 40) are the same used in Ref. |4j for analogous calculations in 
solid p-H 2 : 



Xn{r) 



{1 





cos 



2nn 
L/2-r c 



L/2)]} 



>r 5 r > r r 



r < r r 



(11) 



The cutoff radius r c allows to remove divergences for r — > and therefore to avoid numer- 
ical instabilities in the optimization procedure. If such cutoff is small enough, it does not 
influence the value of the energy. In our calculations it has been taken equal to lA. The 
parameters b, C and {a n } have been determined by following a reweighting scheme, which 
provides to alternate between VMC calculations of (Ti) and minimization of a combination 
of the variance of the expectation value of the Hamiltonian, and of the expectation value 
itself on a set of sampled configurations. 

For each density studied, DMC calculations have been performed projecting from the 
optimized trial function using different imaginary time-steps and populations of walkers, 
in order to check that the bias due to the use of a finite time-step and a limited number 
of walkers were both negligible compared to statistical errors affecting the results of the 
simulations. As well as the ground state energy, the configuration generated in the QMC 
simulations asymptotically sampling from ^t4>o can be used to evaluate the mixed estimator 

{ ~ (*r|0o> " N w ^ * T (Ri) { ] 

of any operator O of interest. The mixed estimator coincides with the ground-state ex- 
pectation value (0o I O\<fio) only if O commutes with the Hamiltonian, otherwise one can 
obtain an estimate of the ground-state expectation by computing the extrapolated estima- 
tor (C) cxt = 2(C) mix - {^ t \0\® t )/{^t\®t) The bias affecting such estimate is of the second 
order in a, where 0o = + Q:^/ a . 

We performed simulations at different densities (0.0194A 3 < p < 0.0430A 3 ) for two 
different lattices, face centered cubic (fee) and hexagonal close packed (hep). Simulation 
cell were set to accommodate 3x3x3 elementary cubic cells for fee, with a total of iV=108 
lattice sites, and 5x3x3 elementary cells for a total of iV=180 lattice sites for hep. Periodic 
boundary conditions are also imposed in order to reduce finite-size effects. The molecule- 
molecule interaction is truncated at the edge of a sphere of radius equal to L/2, where L is 
the length of the shortest side of the cell. In order to avoid discontinuities in the Green's 
function used to sample the walkers, the potential was also shifted of a quantity v(L/2). The 
contribution from the potential energy outside the sphere is estimated by integrating the 
potential in the (L/2, +oo] interval, assuming therefore a uniform pair correlation function 
of molecules beyond L/2. A similar procedure has been used for correcting the error due to 
the shift in the sphere of radius L/2. 



i 




-125 



• SG fee 

.120 " SG *<* - 

♦ Buck /cc ft 

a Buck ftcp 1 if 
135 * exp /,< 

// 



1016 0.020 0.024 0.028 0.032 0.036 0.040 0.044 
P (A 3 ) 



FIG. 1: Energy per particle (in K) in solid P-H2 as a function of the density p, computed with 
different interactions. The curves represent the fits to the DMC results with Murnaghan-like curve. 
In the inset the region around the equilibrium density is expanded. The experimental point is taken 
from Ref. 
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RESULTS 

In Fig. [I] we report the results obtained for the equation of state in the density interval 
0.0194A" 3 < p < 0.0430A" 3 , for the SG and the Buck potentials. In the inset we report 
the expanded curves for density near the equilibrium one. The curves report the DMC 
results, which are found to be only slightly lower than the corresponding results obtained 
by optimizing the Jastrow factor in the wavefunction. We found that the difference between 
VMC and DMC results is less than 2.5K in all the range of density studied and depends 
only weakly on the density. This suggests that the trial wavefunction used is already quite 
accurate. 

We fitted the computed DMC results of the energy per particle (Table [I]) by means of a 
mo died Murnaghan curve: 

e(p) =a + bp + cp\ (13) 

The coefficients and the minimum of the curves in the fee and hep crystal are reported 
in Table IIHI The only available experimental estimate of the sublimation heat at OK for 
solid D 2 is -135. 5K per molecule 6| and refers to solid normal D 2 (33% p-D 2 ), whose binding 
energy is expected to be slightly lower than that of o-D 2 . The experimental value can 
be therefore considered only a lower bound for the binding energy of pure o-D 2 . In the 
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p 


SG fee 


SG hep 


Buck fee 


Buck hep 


0.01940 


-97.21(1) 


-97.37(3) 


-101.73(2) 


-101.91(2) 


0.02150 


-108.23(2) 


-108.48(4) 


-113.42(2) 


-113.69(1) 


0.02509 


-123.70(2) 


-124.09(2) 


-129.99(1) 


-130.38(2) 


0.02609 


-126.89(2) 


-127.33(2) 


-133.47(1) 


-133.90(1) 


0.02900 


-132.55(2) 


-133.08(2) 


-139.81(2) 


-140.33(1) 


0.03170 


-132.05(2) 


-132.67(2) 


-139.73(1) 


-140.35(1) 


0.03400 


-126.49(3) 


-127.14(1) 


-134.30(1) 


-135.03(2) 


0.03700 


-111.11(3) 


-111.97(2) 


-118.77(1) 


-119.62(2) 


0.04000 


-85.59(3) 


-86.50(3) 


-92.58(1) 


-93.50(2) 


0.04300 


-48.89(4) 


-49.92(3) 


-54.63(3) 


-55.71(2) 



TABLE I: Energy per particle (in K) in solid 0-D2 computed with different interaction for fee and 
hep crystalline structure. 



p 


108 (SG) 


256 (SG) 


108 (Buck) 


256 (Buck) 


0.02609 


-126.89(2) 


-126.77(1) 


-133.47(1) 


-133.35(1) 


0.04300 


-48.89(4) 


-48.52(2) 


-54.63(3) 


-54.33(3) 



TABLE II: Energy per particle (in K) in solid 0-D2 computed in fee crystal using different number 
of molecules interacting via the Buck potential. Densities are in A~ 3 . 





a 


b 


c 


7 


Po 


eo 


SG fee 


80.4313 


-9932.81 


1.81514-10 7 


3.50150 


0.03009(1) 


-133.12(4) 


SG hep 


79.9612 


-9912.63 


1.86138-10 7 


3.51105 


0.03014(1) 


-133.69(4) 


Buck fee 


85.0026 


-10397.5 


2.0420M0 7 


3.52878 


0.03028(1) 


-140.63(4) 


Buck hep 


83.5366 


-10313.4 


2.16255-10 7 


3.55035 


0.03033(1) 


-141.14(4) 



TABLE III: Coefficients of the equation of state (Eq. [TBI) fitted to DMC simulation results and 
computed equilibrium energy and equlibrium density. The equilibrium energies eo are given in 
Kelvin; the equilibrium densities po are given in A~ 3 . Units of the parameters a, b and c are 
expressed accordantly. 
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case of molecular hydrogen, the difference in energy between normal and pure P-H2 is less 
than 4%0, 6]. The computed DMC values for 0-D2 at equilibrium density of 0-D2 using 
the Buck potential (E/N = -140.63K and E/N = -141. 14K for the fee and hep crystal 
respectively) are found to be about 5K below the experimental result. Such result confirms 
that employing an effective potential which completely neglects tree-body interactions, like 
the Buck one, leads to an overestimate of the binding energy of the system, as found in the 
case of 4 Hefl and p-H 2 fl. On the other hand the SG potential gives results about 2K above 
the experimental lower bound. This discrepancy is due to the effective triple-dipole term 
proportional to 1/r 9 , which is repulsive, and which therefore tends to increase the energy 
per particle. If we suppose that this effective term tends to overestimate the equilibrium 
energy as in the case of p-H 2 , we can estimate the difference between the value found for 
the GS potential and the right value to be less than 2%. 

In both the cases the hep lattice turns out to be always stable with respect to the fee one. 
The equilibrium density computed by means of the SG potential has been found in better 
agreement with experimental finding po = 0.03009A 3 of Ref. s|. 

In Tab. [IT] we also reported some results of simulations performed in the fee crystal using 
256 molecules arranged on 4x4x4 elementary cells. The bias affecting the estimates of the 
energy due to the finite size of the simulation box have been found negligible in comparison 
with the energetic differences provided by the choice of different potentials. The energy gap 
between fee and hep lattice seems to slightly increase if the fee lattice with 256 sites is taken 
into account. 

The pressure as a function of the density was then obtained from the fitted curves by 
means of the following expression: 



P{p) = P 



dp 



(14) 



The computed curves for the pressure are reported in Fig. [21 and compared with the 
experimental results of Ref. j& for pure o-D 2 . At low densities the curves obtained the two 
potentials tend to overlap and provide only a slight underestimation of the experimental 
pressure. For densities p > 0.037A~ 3 the curve given by the Buck potential tends clearly 
to overestimate the pressure, while the SG potential seems to be capable to keep a much 
better agreement with experimental data, although the density dependence of the pressure 
indicates that the pressure would tend to be overestimated at higher densities. It can be 
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FIG. 2: Pressure as a function of the density p in solid P-H2 computed with different interactions 
Curves for fee and hep are almost indistinguishable. The experimental points are taken from Ref. 

expected that at high densities the effects of many body terms in the potential become more 
important, but the the triple-dipole effective term in the SG potential seems to be capable 
to take into account most of the many-body effects on the pressure curve. 

The energy gap between the Buck and SG potential can be therefore considered as a first 
estimate of the contribution of three body terms to the potential, That contribution, ranging 
from 3 to 7K, represents about the 4% of the absolute vaule of the two-body potential energy 
in all the range of studied density. In the case of solid p-H 2 , indeed, the contribution of 
the three body term, although comparable, represents from 4% up to the 12% of the two- 
body potential energy, depending on the density, leading to much stronger effects on the 
P-V curve, which have to be kept into account by adding an explicit Axilrod-Teller three- 
body term to the potential [4|. The smaller relative influence of the many-body terms on 
the equation of state of o-D 2 is to attribute to its heavier mass, which tends to keep the 
molecules more localized and therefore to increase the two-body potential energy. 

The many-body effects in o-D 2 are stronger than of 4 He, although they have the same 
mass: the energetic contribution of the three-body potential in 4 He was estimated to be 
less than 1.2% of the absolute value of the two-body potential energy at densities from 
0.02934A' 3 to 0.03527A- 3 . 

Once the equation of state e(p) is known, it is possible to calculate the isothermal com- 
pressibility, defined as: 

1 [ An~\ 

(15) 



K( P ) = - 



dp 
dP 
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FIG. 3: Compressibility in molecular P-H2 computed using the SG and the Buck interactions. 
Curves for fee and hep are almost indistinguishable 




r(A) 



r(A) 



FIG. 4: Pair distribution function in solid p-H2Computed using the SG and the Buck potentials at 
density p=0.04300A -3 both in the fee and in the hep crystal. 

In Fig. [3] the results for the compressibility yelded by the two potentials are compared 
with the corresponding experimental value obtained by fitting the pressure curve of Ref. 
Isl. Both potentials provide a rather good agreement with experimental data except at the 
lowest densities, as expected due the bigger relative discrepancy between experimental and 
calculated values of the pressure that arises when the pressure tends to vanish. In order to 
investigate the effects of the choice of different models of intermolecular potentials on the 
local structure of the o-D 2 solid we computed the pair distribution function: 



9(r) 



i+3 



(16) 



Results for the computation of the function g(r) are reported in Fig. HI For all the 
range of density studied the curves obtained employing the SG and the Buck potential are 
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almost indistinguishable from each other, giving an indication that the effects of the choice 
of different models of interaction on the structural properties of the crystal of o-D 2 are 
negligible, as already found also for solid p-H 2 . 

CONCLUSIONS 

We present the results of a systematic comparison of ground state properties of solid o-D 2 
computed with the two most used interactions, the Buck and the Silvera-Goldman potential. 
Only the last one is capable to give a satisfactory description of the pressure-volume curve. 
We attribute this result to the presence in the potential of a repulsive extra-term which 
treats in an effective way the possible three body interactions acting among the molecules. 

No evident structural differences have been found by changing interaction model. 
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